#coding:utf-8

from math import sqrt

n=5
xi = [0,1,2,3,4]
yi = [3,4,1,5,2]
zi = [2,3,3,1,5]
ri = [6, 7, 2, 3, 4]

def soma_v(v1,v2):
	res = [0]*len(v1)
	for i in range(0,len(v1)):
		res[i] = v1[i]+v2[i]
	return res
	
def interno_v(v1,v2):
	res = 0
	for i in range(0,len(v1)):
		res += v1[i]*v2[i]
		
	return res	
	
def escalar_v(v1,k):
	res = [0]*len(v1)
	for i in range(0,len(v1)):
		res[i] = v1[i]*k
		
	return res	
	
def norma(v):
	res = 0
	for i in range(0,len(v)):
		res+=pow(v[i],2)
	res = sqrt(res)
	return res
	
def Ei(i,k,v):
	ind_i = 3*i
	ind_k = 3*k
	
	return  pow( (pow( (v[ind_k]-v[ind_i]),2 ) + pow( (v[ind_k+1]-v[ind_i+1]),2 ) + pow( (v[ind_k+2]-v[ind_i+2]),2 ) - pow( (ri[k]+ri[i]),2 ) ),2 )

def der_Ei_k(i,k,v):
	ind_i = 3*i	
	ind_k = 3*k
	
	dx = 4*(pow( (v[ind_k]-v[ind_i]),2 ) + pow( (v[ind_k+1]-v[ind_i+1]),2 ) + pow( (v[ind_k+2]-v[ind_i+2]),2 ) - pow( (ri[k]+ri[i]),2 ) )*(v[ind_k]-v[ind_i])
	dy = 4*(pow( (v[ind_k]-v[ind_i]),2 ) + pow( (v[ind_k+1]-v[ind_i+1]),2 ) + pow( (v[ind_k+2]-v[ind_i+2]),2 ) - pow( (ri[k]+ri[i]),2 ) )*(v[ind_k+1]-v[ind_i+1])
	dz = 4*(pow( (v[ind_k]-v[ind_i]),2 ) + pow( (v[ind_k+1]-v[ind_i+1]),2 ) + pow( (v[ind_k+2]-v[ind_i+2]),2 ) - pow( (ri[k]+ri[i]),2 ) )*(v[ind_k+2]-v[ind_i+2])
	
	return [ dx,dy,dz ]
	
def Fi(i,v):
	ind = 3*i
	res = 0
	for k in range(0,n):
		if k!=i:
			res+=Ei(k,i,v)
	
	return res
	
def der_Fi_k(i,k,v):
	ind_i=3*i
	ind_k=3*k
	
	res = [0]*3
	
	if i==k:
		for j in range(0,n):
			if j!=i:
				soma_v(res,der_Ei_k(j,i,v))
	else:
		dx = -4*(pow( (v[ind_k]-v[ind_i]),2 ) + pow( (v[ind_k+1]-v[ind_i+1]),2 ) + pow( (v[ind_k+2]-v[ind_i+2]),2 ) - pow( (ri[k]+ri[i]),2 ) )*(v[ind_i]-v[ind_k])
		dy = -4*(pow( (v[ind_k]-v[ind_i]),2 ) + pow( (v[ind_k+1]-v[ind_i+1]),2 ) + pow( (v[ind_k+2]-v[ind_i+2]),2 ) - pow( (ri[k]+ri[i]),2 ) )*(v[ind_i+1]-v[ind_k+1])
		dz = -4*(pow( (v[ind_k]-v[ind_i]),2 ) + pow( (v[ind_k+1]-v[ind_i+1]),2 ) + pow( (v[ind_k+2]-v[ind_i+2]),2 ) - pow( (ri[k]+ri[i]),2 ) )*(v[ind_i+2]-v[ind_k+2])
		
		res[0] = dx
		res[1] = dy
		res[2] = dz
		
	return res	
	
def F(v):
	res = 0
	
	max = [-100,-100,-100]
	min = [100,100,100]		
	for i in range(0,n):
		res+=Fi(i,v)
		ind = 3*i
		
		if v[ind]+ri[i] > max[0]:
			max[0] = v[ind]+ri[i]
		if v[ind+1]+ri[i] > max[1]:
			max[1] = v[ind+1]+ri[i]
		if v[ind+2]+ri[i] > max[2]:
			max[2] = v[ind+2]+ri[i]			
			
		if v[ind]-ri[i] < min[0]:
			min[0] = v[ind]-ri[i]
		if v[ind+1]-ri[i] < min[1]:
			min[1] = v[ind+1]-ri[i]
		if v[ind+2]-ri[i] < min[2]:
			min[2] = v[ind+2]-ri[i]
		
	res += pow( ( max[0]-min[0] ),2) + pow( (max[1]-min[1]),2 ) + pow( (max[2]-min[2]),2 )		
		
	return res
	
def der_F_k(k,v):
	res = [0]*3
	for i in range(0,n):
		dfik = der_Fi_k(i,k,v)
		res[0]+=dfik[0]
		res[1]+=dfik[1]
		res[2]+=dfik[2]
		
	return res	
	
def grad_F(v):
	res = [0]*3*n
	
	max = [-100,-100,-100]
	min = [100,100,100]		
	xk_max = None
	xk_min = None
	yk_max = None
	yk_min = None
	zk_max = None
	zk_min = None
	
	for k in range(0,n):
		ind = 3*k
		dfk = der_F_k(k,v)
		res[ind] = dfk[0]
		res[ind+1] = dfk[1]
		res[ind+2] = dfk[2]
		
		if v[ind]+ri[k] > max[0]:
			max[0] = v[ind]+ri[k]
			xk_max = k
		if v[ind+1]+ri[k] > max[1]:
			max[1] = v[ind+1]+ri[k]
			yk_max = k
		if v[ind+2]+ri[k] > max[2]:
			max[2] = v[ind+2]+ri[k]			
			zk_max = k
			
		if v[ind]-ri[k] < min[0]:
			min[0] = v[ind]-ri[k]
			xk_min = k
		if v[ind+1]-ri[k] < min[1]:
			min[1] = v[ind+1]-ri[k]
			yk_min = k
		if v[ind+2]-ri[k] < min[2]:
			min[2] = v[ind+2]-ri[k]		
			zk_min = k	
		
		if xk_min is not None:
			res[xk_min]+= -2*(max[0]-min[0])
		if yk_min is not None:
			res[yk_min]+= -2*(max[1]-min[1])
		if zk_min is not None:
			res[zk_min]+= -2*(max[2]-min[2])	
			
		if xk_max is not None:
			res[xk_max]+= 2*(max[0]-min[0])
		if yk_max is not None:
			res[yk_max]+= 2*(max[1]-min[1])
		if zk_max is not None:
			res[zk_max]+= 2*(max[2]-min[2])
		
	return res
	
def v_from_vec():
	v = [0]*3*n
	for i in range(0,n):
		ind = 3*i
		v[ind] = xi[i]
		v[ind+1] = yi[i]
		v[ind+2] = zi[i]
		
	return v

def busca_linear(xk,dk):
	sigma = 0.8
	lbda = 1.0
	lim = 10
	
	while( F( soma_v(xk,escalar_v(dk,lbda)) ) >= F(xk) + sigma*lbda*interno_v(grad_F(xk),dk) ):
		if lim == 0:
			break			
		lbda/=2
		lim-=1
	
	return lbda
		
def main():
	x0 = v_from_vec()
	epson = 1e-2
	it = 0
	
	while(norma(x0)>epson):
		if it==100:
			break
			
		dk = escalar_v(grad_F(x0),-1)
		lbda = busca_linear(x0,dk)
		x0 = soma_v(x0,escalar_v(dk,lbda))
		
		it+=1
	
	for i in range(0,n):
		ind = 3*i
		print "P",i,"= ( ",x0[ind],x0[ind+1],x0[ind+2],ri[i]," )"
		
	return x0
	
#x = main()
x0 = v_from_vec()